#pragma once
#include <zCoreConfig.hpp>
#include "../common.hpp"
#include "Vector3.hpp"
#include "Matrix4x4.hpp"

namespace zzz{
#undef max
#undef min
inline zuint TimeSeed() {
  return time(NULL);
}

inline void RandSeed() {
  srand(TimeSeed());
}


template<typename T>
T UniformRand(const T &min, const T &max) {
  int r=rand();
  double rr=(double)r/(double)RAND_MAX;
  return (max-min)*rr+min;
}

// The Park and Miller "Minimal Standard" LCG generator
// X[n+1] = (X[n]*a)%m
// Adapted from Numerical Recipies in C, 2nd Ed
// Similar to FreeBSD random number generator
class ZCORE_CLASS RandomLCG {
public:
  RandomLCG(const zuint seed = 199125);  ///< Construct with non-zero seed
  RandomLCG(const RandomLCG &rng);    ///< Copy constructor
  zuint Rand();                        ///< Random number in [0,max()]
  zuint Max()  const;                        ///< Maximum random number
  void Seed(const zuint seed = 199125);      ///< Seed the random number generator
  void SeedFromTime();
private:
  mutable zuint idum_;
  const static int a;
  const static int m;
  const static int q;
  const static int r;
};

// Multiple LFSR (Linear Feedback Shift Register) generator
// Larry Smith
// Nigel Stewart, RMIT (nigels@nigels.com)
// Original Pascal implementation from unknown source
// This random number generator uses three LFSRs (Linear Feedback Shift Registers)
// to create long sequences of random bits.  Each bit is equally randomised.
// The origins of this exact algorithm is unknown.  This implementation is
// based on the ctools module rndnums by Larry Smith, used with permission.
// It is not really known how good this method is for random number generation -
// but it has the advantage of generating 32 random bits.  Probably not recommended
// for crypto applications....
// Further information about LFSRs:
// B. Schneier, Applied Cryptography, 2nd Ed, Johyn Wiley & Sons
class ZCORE_CLASS RandomLFSRMix {
public:
  /// Constructor
  RandomLFSRMix(const zuint seed1 = 199125,const zuint seed2 = 90618,const zuint seed3 = 189419);
  RandomLFSRMix(const RandomLFSRMix &rng);  ///< Copy constructor
  zuint Rand();                            ///< 32-bit random number
  zuint Max()  const;                            ///< Maximum random number
  /// Seed the random number generator
  void Seed(const zuint seed1 = 199125,const zuint seed2 = 90618,const zuint seed3 = 189419);
  void SeedFromTime();
  static RandomLFSRMix rng;                    ///< Global LFSRMix RNG for convenience
private:
  mutable zuint regA_, regB_, regC_;
};

// Adaptor for random numbers in integer domain
template<typename T, class R = RandomLFSRMix>
class RandomInteger {
public:
  /// Constructor
  RandomInteger(R &random,const T min,const T max)
  :random_(random),min_(min),mod_(max-min+1) { ZCHECK_LE(min, max); }
  /// Constructor using default RNG
  RandomInteger(const T min,const T max)
  :random_(R::rng),min_(min),mod_(max-min+1) { ZCHECK_LE(min, max); }
  /// Copy constructor
  RandomInteger(const RandomInteger<T, R> &gen) : random_(gen.random_),min_(gen.min_),mod_(gen.mod_){}
  /// Random number in the range [min,max]
  zuint Rand() {
    return min_ + random_.Rand()%mod_;
  }
  void SeedFromTime(){random_.SeedFromTime();}
private:
  R& random_;
  const T min_;
  const T mod_;
};

// Adaptor for random float or double numbers
template<typename T ,class R = RandomLFSRMix>
class RandomReal {
public:
  /// Constructor
  RandomReal(R &random,const T min = 0.0,const T max = 1.0)
    : random_(random),min_(min),max_(max),mult_(max-min){}
  /// Constructor using default RNG
  RandomReal(const T min = 0.0,const T max = 1.0)
    : random_(R::rng),min_(min),max_(max),mult_(max-min){}
  /// Copy constructor
  RandomReal(const RandomReal<T, R> &gen)
    : random_(gen._random),min_(gen.min_),max_(gen.max_),mult_(gen.mult_){}

  /// Random number in the range [min,max]
  T Rand() {
    return Lerp(min_,max_,float(random_.Rand())/numeric_limits<zuint>::max());
  }
  void SeedFromTime(){random_.SeedFromTime();}
private:
  R& random_;
  const T min_;
  const T max_;
  const T mult_;
};

//normal distribution (0, 1)
template<typename T, class R = RandomLFSRMix>
class RandomNormal {
public:
  /// Constructor
  RandomNormal(R &random,const T mean = 0, const T var = 1) : z_(random, 0.0, 1.0),mean_(mean),var_(var),takenX2_(true){}
  /// Constructor using default RNG
  RandomNormal(const T mean = 0, const T var = 1) : z_(R::rng, 0.0, 1.0),mean_(mean),var_(var),takenX2_(true){}
  /// Copy constructor
  RandomNormal(const RandomNormal<T, R> &gen) : z_(gen.z_.random_),mean_(gen.mean_),var_(gen.var_),takenX2_(true){}
  /// Random point on sphere surface
  T Rand() {
    //Box-Muller method
    if (takenX2_==false) {
      takenX2_=true;
      return lastX2_;
    }
    T R1=z_.Rand(),R2=z_.Rand();
    lastX2_=Sqrt(-2*log(R1))*sin(C_2PI*R2)*var_+mean_;
    takenX2_=false;
    return Sqrt(-2*log(R1))*cos(C_2PI*R2)*var_+mean_;
  }
  void SeedFromTime(){z_.SeedFromTime();}
private:
  RandomReal<T, R> z_;
  const T mean_, var_;
  bool takenX2_;
  T lastX2_;
};

// Random points on the unit sphere
// Nigel Stewart, RMIT (nigels@nigels.com)
// http://www.math.niu.edu/~rusin/known-math/96/sph.rand
// The trig method.  This method works only in 3-space, but it is
// very fast.  It depends on the slightly counterintuitive fact
// that each of the three coordinates of a uniformly
// distributed point on S^2 is uniformly distributed on [-1, 1] (but
// the three are not independent, obviously).  Therefore, it
// suffices to choose one axis (Z, say) and generate a uniformly
// distributed value on that axis.  This constrains the chosen point
// to lie on a circle parallel to the X-Y plane, and the obvious
// trig method may be used to obtain the remaining coordinates.
// 
// Choose z uniformly distributed in [-1, 1].
// Choose t uniformly distributed on [0, 2*pi).
// Let r = sqrt(1-z^2).
// Let x = r * cos(t).
// Let y = r * sin(t).
template<typename T, class R = RandomLFSRMix>
class RandomSphere {
public:
  /// Constructor
  RandomSphere(R &random,const T radius = 1.0) : z_(random, -1.0, 1.0),t_(random, 0,C_2PI),radius_(radius){}
  /// Constructor using default RNG
  RandomSphere(const T radius = 1.0) : z_(R::rng, -1.0, 1.0),t_(R::rng, 0.0,C_2PI),radius_(radius){}
  /// Copy constructor
  RandomSphere(const RandomSphere<R> &gen) : z_(gen.z_.random_),t_(gen.t_.random_),radius_(gen.radius_){}
  /// Random point on sphere surface
  Vector<3, T> Rand()
  {
    const T z = z_.Rand();
    const T t = t_.Rand();
    const T r = radius_*Sqrt<T>(1.0-z*z);
    return Vector<3, T> (r*cos(t), r*sin(t),radius_*z);
  }
  void SeedFromTime(){z_.SeedFromTime();t_.SeedFromTime();}
private:
  RandomReal<T, R> z_,t_;
  const T radius_;
};

// Random orientation matricies in 3D
template<typename T, class R = RandomLFSRMix>
class RandomOrientation {
public:
  /// Constructor
  RandomOrientation(R &random):s_(random),a_(random, 0.0,C_2PI){}
  /// Constructor using default RNG
  RandomOrientation():s_(R::rng),a_(R::rng, 0.0,C_2PI){}
  /// Copy constructor
  RandomOrientation(const RandomOrientation<R> &gen):s_(gen.s_),a_(gen.a_){}
  /// Random orientation matrix
  Matrix<4, 4, T> Rand() {
    // A is a vector distributed
    // randomly on the unit sphere
    const Vector<3, T> a(s_.Rand());
    // B is orthogonal to A
    Vector<3, T> b;
    switch (a.MaxPos()) {
    case 0: b = a.Cross(Vector<3, T>(0, 0, 1)); break; // a in x direction
    case 1: b = a.Cross(Vector<3, T>(0, 0, 1)); break; // a in y direction
    case 2: b = a.Cross(Vector<3, T>(1, 0, 0)); break; // a in z direction
    }
    // C is orthogonal to both A and B
    Vector<3, T> c(b.Cross(a));
    // TODO - Possible to avoid this step?
    b.Normalize();
    c.Normalize();
    // Choose an angle in the AB plane
    const T angle = a_.Rand();
    const T sina  = sin(angle);
    const T cosa  = cos(angle);

    return Matrix<4, 4, T>(a, b*sina+c*cosa, b*cosa-c*sina);
  }
  void SeedFromTime(){s_.SeedFromTime();a_.SeedFromTime();}
private:
  RandomSphere<T, R> s_;
  RandomReal<T, R> a_;
};

//random sampling on unit hypersphere surface
template<int N, typename T, class R = RandomLFSRMix>
class RandomHyperSphere {
public:
  /// Constructor
  RandomHyperSphere(R &random, const T radius = 1.0) : theta_(random, 0.0,C_2PI), radius_(radius){}
  /// Constructor using default RNG
  RandomHyperSphere(const T radius = 1.0) : theta_(RandomLFSRMix::rng, 0.0,C_2PI),radius_(radius){}
  /// Copy constructor
  RandomHyperSphere(const RandomHyperSphere<N, T, R> &gen) : theta_(gen.theta_),radius_(gen.radius_){}
  /// Random point on sphere surface
  Vector<N, T> Rand() {
    Vector<N-1, T> thetas;
    for (int i=0; i<N-1; i++) thetas[i]=theta_.Rand();
    Vector<N, T> res;
    res[0]=1;
    for (int i=1; i<N-1; i++)
      res[i]=res[i-1]*sin(thetas[i-1]);
    res[N-1]=res[N-2];
    for (int i=0; i<N-1; i++)
      res[i]*=cos(thetas[i]);
    res[N-1]*=sin(thetas[N-2]);
    if (radius_==1) return res;
    else return res*radius_;
  }
  Vector<N, T> RandOnDim(int n) {
    assert(n<=N);
    Vector<N-1, T> thetas;
    for (int i=0; i<n-1; i++) thetas[i]=theta_.Rand();
    Vector<N, T> res;
    res[0]=1;
    for (int i=1; i<n-1; i++)
      res[i]=res[i-1]*sin(thetas[i-1]);
    res[n-1]=res[n-2];
    for (int i=0; i<n-1; i++)
      res[i]*=cos(thetas[i]);
    res[n-1]*=sin(thetas[n-2]);
    for (int i=n; i<N; i++)
      res[i]=0;
    if (radius_==1) return res;
    else return res*radius_;
  }

  void SeedFromTime(){theta_.SeedFromTime();}
private:
  RandomReal<T, R> theta_;
  const T radius_;
};

//random sampling on unit hypersphere surface
template<int N, typename T, class R = RandomLFSRMix>
class RandomHyperSphere2 {
public:
  /// Constructor
  RandomHyperSphere2(R &random, const T radius = 1.0) : z_(random, 0.0, 1.0), radius_(radius){}
  /// Constructor using default RNG
  RandomHyperSphere2(const T radius = 1.0) : z_(RandomLFSRMix::rng, 0.0, 1.0),radius_(radius){}
  /// Copy constructor
  RandomHyperSphere2(const RandomHyperSphere<N, T, R> &gen) : z_(gen.z_),radius_(gen.radius_){}
  /// Random point on sphere surface
  Vector<N, T> Rand() {
    Vector<N, T> z;
    for (int i=0; i<N; i++) z[i]=z_.Rand();
    T r=FastDot(z, z);
    r=Sqrt(r);
    Vector<N, T> ret;
    for (int i=0; i<N; i++) ret[i]=z[i]/r;
    T len=ret.Len();
    if (len==0) return Rand();
    ret/=len;
    if (radius_==1) return ret;
    return ret*radius_;
  }
  Vector<N, T> RandOnDim(int n) {
    assert(n<=N);
    Vector<N, T> z;
    for (int i=0; i<n; i++) z[i]=z_.Rand();
    T r=0;
    for (int i=0; i<n; i++) r+=z[i]*z[i];
    r=Sqrt(r);
    Vector<N, T> ret;
    for (int i=0; i<n; i++) ret[i]=z[i]/r;
    for (int i=n; i<N; i++) ret[i]=0;
    T len=ret.Len();
    if (len==0) return RandOnDim(n);
    ret/=len;
    if (radius_==1) return ret;
    else return ret*radius_;
  }
  void SeedFromTime(){z_.SeedFromTime();}
private:
  RandomNormal<T, R> z_;
  const T radius_;
};

}
